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SUMMARY 


An iterative transformation procedure suggested by H. Wielandt for 
numerical solution of flutter and similar characteristic -value problems 
is presented. Application of this “'procedure to ordinary natural- 
vibration problems and to flutter problems is shown by numerical 
examples. Comparisons of computed results with experimental values and 
with results obtained by other methods of analysis are made. 


INTRODUCTION 


Existing methods of flutter analysis include the representative - 
section method, generalized-coordinate methods, matrix methods, and 
operational methods. The present paper presents an iteration procedure 
for analysis of flutter and similar characteristic -value problems. 

For ordinary natural-vibration problems, iterative procedures of 
the Stodola type (references 1 and 2) are suitable for finding the funda- 
mental and higher-order natural modes and frequencies. The higher-order 
solutions are obtained through use, of the orthogonality relations that 
exist among the natural modes. 

Orthogonality relations analogous to those that exist in ordinary 
vibration problems can be found in the flutter problem only by intro- 
duction of the so-called "adjoint" problem. (This additional step is 
unnecessary in ordinary vibration problems by virtue of the fact that 
they are "self-adjoint.") Wielandt has suggested an iterative transfor- 
mation procedure (reference 3) which is well -suited to the flutter 
problem in that it avoids the need of orthogonality relations and hence 
does not require consideration of the adjoint problem. The iterative 
transformation procedure can also be applied to ordinary natural- 
vibration problems with less labor than is generally required in the 
extended Stodola procedure. 
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Because both the original and translated works of Wielandt are 
difficult to follow, an explanation of the idea of the iterative trans- 
formation procedure is given in the present paper and the application 
of the procedure to ordinary natural -vibration problems and to flexure - 
torsion flutter problems is shown in numerical examples. Comparisons 
of computed results with experimental values and with results obtained 
by other methods of analysis are also made . 


SYMBOLS 


El 

GJ 

x 

y 

0 


w 


P 0 


V S 0 


flexural stiffness 

torsional stiffness 

spanwise coordinate with origin at root of wing 

complex representation of amplitudes and phases of translation 
of elastic axis in harmonic motion 

/ 

complex representation of amplitudes and phases of rotation 
about elastic axis in harmonic motion 

coupled mode (y,0) 

complex coefficients of y which, when multiplied by y, give 
complex representation of amplitudes and phases of aero- 
dynamic and inertia forces associated with translational 
component of harmonic motion 

complex coefficients of 0 which, when multiplied by 0, give 
complex representation of amplitudes and phases of aerodynamic 
and inertia forces associated with rotational component of 
harmonic motion 

complex coefficients of y which, when multiplied by y, give 
complex representation of amplitudes and phases of aerodynamic 
and inertia torques associated with translational component 
of harmonic motion 

complex coefficients of 0 which, when multiplied by 0, give 
complex representation of amplitudes and phases of aerodynamic 
and inertia torques associated with rotational component of 
harmonic motion 

structural -damping coefficients associated with flexure and 
torsion, respectively (see appendix B) 
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g a coefficient of artificial damping (may be either positive or 

negative) 

k reduced frequency (bcu/v) 

0) frequency of harmonic motion 

C characteristic value 

b length of semichord of wing 

L length of cantilever wing from root to tip 

P mass ratio (7/jtpb^) 

v • velocity of air relative to wing 

7 distributed mass of wing per unit length of span 

p mass, density of air 

a distance between midchord axis and elastic axis in terms of 

local semichord, positive when elastic axis is behind mid- 
chord axis 

u distance between elastic axis and gravity axis of distributed 

mass of wing in terms of local semichord, positive when 
gravity axis is behind elastic axis 

r radius of gyration of distributed mass of wing about elastic 

axis in terms of local semichord 

F,G transcendental functions of k (see reference 4) 

t time 



F Tnn eigenvalue factor Up - lj 

R ratio of complex constants 

X length; in numerical solutions, distance between specific 

adjacent stations of wing 


P 


applied force 
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q applied torque 

V shear 

M bending moment 

a curvature 

3 slope of elastic axis 

T twisting moment 

9 angle of twist 

Subscripts: 

' fi 

true modes or eigenvalues 
a2,a3,a4, ... transformed modes 

b intermediate derived mode 

A,B,C,... stations 

R real 

I imaginary 

o reference value 

bl,ba2,ba3, . . . sweeping functions 
Superscripts: 

(1),(2),(3),.. . cycles of iteration 

A bar over a symbol indicates a concentrated quantity instead of a 
distributed quantity. 

A prime is used to denote division by a £ . 
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ITERATIVE TRANSFORMATION METHOD OF SOLUTION 
General Features of Method 


The principle of the iterative transformation procedure is similar 
in form to that of the standard iteration procedure for solving 
characteristic -value problems. Both procedures require the determina- 
tion of the solutions in the order of the magnitudes of the eigenvalue s, 
beginning with the fundamental. Both procedures require assumptions of 
modes, integrations which generally must be done numerically, and 
sweeping operations for higher -order -mode determinations. The distin- 
guishing features of the iterative transformation procedure occur in 
the determination of solutions higher than the fundamental and are as 
follows: (l) The immediate aim is to determine not the true nth mode, 

as in the standard iteration procedure, but a particular linear combina- 
tion composed of all modes from the fundamental to the nth. This linear 
combination is referred to as the transformed nth mode. The transformed 
nth mode can be made to have nodal (zero) points at specified stations 
of the wingj such a feature is highly desirable in numerical work. 

( 2 ) The sweeping operations, which consist of subtractions of lower- 
order-mode shapes from the function obtained by integrating the assumed 
mode, do not employ the orthogonality relations as in the standard 
iteration method but make use of forcing functions that, in numerical 
work, greatly simplify the sweeping operations and increase the over- 
all accuracy of the results by making the sweeping operations more 
consistent with the rest of the process. (3) Although the true nth 
eigenvalue is determined directly in the iterative transformation pro- 
cedure, the true nth mode must be computed from quantities within the 
iteration cycle after the transformed nth mode is found. 


Outline of Steps in the Procedure 


The equation of equilibrium of a cantilever beam vibrating harmoni- 
cally in pure flexure is 


d^ _ T d^y 2 
— - El — | = 7u> y 


-EI i 
dx 1- dx" 

or, after integration with proper attention to boundary conditions, 


( 1 ) 




( 2 ) 
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The solutions of this integral equation are the true natural modes 
(eigenfunctions) y^ , y^, y^, . . . and the corresponding natural 

frequencies (eigenvalues) o^, o^, CD 3 , .... For convenience in 

subsequent discussion, the true modes axe assumed to be normalized to 
unity at some position (station A) along the beam. 

The first mode and frequency are assumed to have been previously 
determined by the Stodola process. The iterative transformation pro- 
cedure becomes applicable in the determination of the second mode and 
frequency. As mentioned previously, the immediate aim in the iteration 
for the second mode. is the determination of a linear combination of first 
and second modes which is called the transformed second mode. The linear 
combination y 2 - y^ which has zero ordinate at station A is chosen and 

defined as the transformed second mode to be determined. The iteration 
for determination of this transformed second mode may be described as 
follows : 

(!) A plausible shape y^ for the transformed second mode is 

assumed. This shape must have zero ordinate at station A and should 
satisfy the boundary conditions as closely as possible. 

(2) The displacement 






resulting from the inertia load 70 ^ y a 2 corresponding to the assumed 
shape y^ 2 ^ • vibrating harmonically at frequency is calculated. 

This calculation must usually be done numerically with the square of the 
frequency being carried along as an undetermined factor. 

(3) A first-mode shape (previously determined) is subtracted 
(swept out) from the calculated displacement y^ in an amount such 

that the resulting displacement is zero at station A. Thus the 
resulting displacement is 



( 2 ) 


(4) The resulting displacement y 

' 1 \ 3*2 

di splacement y^ 


is compared with the assumed 


When the computations are numerical y the 
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ratios y a2 / y a 2 8X6 com P are <i at all the stations. If the assumed 

displacement is exactly equal to the transformed second mode, the ratios 
are equal to each other. These ratios contain the single unknown cup, 
and the second frequency is that value of cup which makes the ratios 
unity . 

(5) If the ratios y^ / ^ a p^ from the first cycle of iteration 

outlined in the four preceding steps are not reasonably the same at nl 1 
stations, the process must be repeated until the ratios become reasonably 
the same. Each new cycle starts with the resultant displacement of each 
preceding cycle. The convergence of this process to the second frequency 
and the transformed second mode is proved in appendix A. 

The transformed third mode and the third frequency are computed in 
the following manner. The transformed third mode is defined as that 
combination of the first three natural modes which has a zero ordinate 
at the same station that was used in the transformed second mode (sta- 
tion A) and also a zero ordinate at some other station, station B. Thus 
the transformed third mode is defined as 


y 3 ' y l " 



The 'iteration is as follows: 

(1) A plausible shape y a ^ for the transformed third mode is 

assumed. This shape must have zero ordinates at stations A and B and 
should satisfy the boundary conditions as closely as possible. 

(2) The displacement 



2 ( 1 ),, .4 
70)3 y a3 '(dx) 


is calculated, with the square of the frequency carried along as 

an undetermined factor. ' 

(3) The first of two sweeping operations, in which a first-mode 
shape is swept from the displacement y b so as to, make the resulting 
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displacement at station A zero, is performed. This operation gives the 
displacement 


y*, - Ur- 


7i/a 


(4) The second sweeping operation, in which a transformed-second- 
mode shape (previously determined) is swept from the resulting displace- 
ment of step (3) so that the new resulting displacement is zero at 
station B as well as at station A, is performed. (This second sweeping 
operation cannot disturb the zero condition at station A established in 
step (3) because the second sweeping function (the transformed second 
mode) is identically zero at station A.) Thus, the final resulting dis- 
placement is 




^b - 




a2 


(5) Comparisons of the ratios y 


(2) /y (1) 
a3 / y a3 


at all stations are made. 


and, if they axe not reasonably the same, additional cycles of iteration 
are carried out until the ratios become reasonably the same. The third 
frequency is then computed from the ratios as explained previously. 
Convergence of this process to the third frequency and the transformed 
third mode is proved in appendix A. 


Frequencies and transformed modes higher than the third may be 
computed by extensions of the process just described. 


Physical Interpretation of the Procedure 

A physical interpretation of the iterative transformation procedure 
can be given. With regard to the transformed second mode, for example, 
the following question may be asked: Under what conditions can the beam 

vibrate in the transformed-second-mode shape at the second natural 
frequency? Vibration in shape y a 2 = yg - yq at frequency implies 

an inertia loading /o>g 2 (yg - yq). But if this load is substituted in 

place of 70 )g 2 y in the right-hand side of equation (2), the result 
after integration will not be y^ - y^ but 
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The inertia and forcing loads are illustrated in figure 1. The inertia 
load acting alone produces a displacement (equation (3)) generally dif- 
ferent from zero at station A. The forcing load produces the displacement 



This displacement (equal to the sweeping function) has the shape of the 
previously determined first mode and is equal and opposite at station A 
to the displacement due to the inertia load) that is, by virtue of the 
previously assigned normalizations at station A, 

($ ^i)* < 6 > 

Thus the displacement due to the forcing load is completely determined 
when the displacement due to the inertia load is known. The gist of the 
foregoing analysis is that vibration in the transformed-second-mode 
shape is the response of the beam to an oscillatory forcing load of the 
first -mode shape and of frequency equal to the second natural frequency, 
superimposed on a free vibration of the beam in the second natural mode. 

Similar physical interpretations of the iterative transformation 
process for modes higher than the second can be made. 


Application of the Procedure in Ordinary 

Coupled Natural -Vibration Problems 

The procedure that has been outlined in a preceding section for 
pure flexure can easily be extended to systems capable of simultaneous 
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flexural and torsional displacements. Airplane wings belong to the 
latter class of systems. The only distinguishing element in coupled 
flexural -torsional vibration problems is that each natural mode contains 
two components, the flexure and the torsion. These components must 
always appear together in a fixed relation to each other. The two com- 
ponents must be computed together and must be used together. 

Each coupled mode is a solution of the simultaneous differential 
equations 


o El ~p = 7<» 2 (y + bu 0 ) 
dx c dx^ 


( 7 ) 


" elx GJ if = ( bu y + b 2 r 2 0 ) (8) 

•Equations ( 7 ) and ( 8 ), after integration, become (for a cantilever beam) 



The solution of the integral equations (9) and (10) for the coupled trans- 
formed second mode by the iterative transformation procedure is outlined 
diagrammatically in figure 2. The flexural component of the displace- 
ment for a particular step is illustrated in the left- hand side of the 
figure and the torsional component is illustrated at the same level in 
the right-hand side. 


In the first step, an approximation to a linear combination of the 
true first and second coupled modes is assumed. The particular linear 
combination having zero flexural displacement at the tip station (sta- 
tion A) is chosen. (For greatest numerical accuracy, this nodal point 
should be chosen in the component and at the station where the first 


coupled mode has its maximum numerical value.) The symbols 

,( 1 ) 

and 0 a2 are used to designate the flexural and torsional 


y (l) 

J a2 

components 


of this assumed displacement, respectively. In general, the magnitude 
of the torsional component relative to the flexural component is 
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difficult to estimate; the most expedient thing to do is take one of 
the components equal to zero. 


The second step is the computation (by numerical integration) of 
the two components of the displacement due to the inertia forces 

7 <» 2 2 (ya 2 + bu 0a2l) 811(1 inertia torques yu ^ 2 (buy^ + b 2 r 2 0 a2 ) that are 
associated with the assumed displacement. The result is termed the 
intermediate derived mode, and the symbols y^ 1 ^ and 0^ are used 

D r D 

to designate its two components. 


third step is the determination of a sweeping function having 
the shape of the first coupled mode (previously determined) and a magni- 
tude such that the sum of the intermediate derived mode and the sweeping 
function equals zero in the flexural component at station A. In 
algebraic terms, the first-mode sweeping function is given by 



( 11 )' 


( 12 ) 


The fourth step is the addition of the intermediate derived mode 
and the first-mode sweeping function to give the derived transformed 
second mode. Thus the two components of the derived transformed second 
mode are 



( 13 ) 



( 14 ) 


The calculation of the ratios 
stations completes the first cycle 


and 

of iteration. 



at all 
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Additional cycles are carried out until the ratios at all stations 
in both the flexural and torsional components have values that are 
reasonably the same. The true second natural frequency of the coupled 
system is then obtained as described previously. 

Steady vibration of an airplane wing at zero airspeed is sin example 
of coupled natural vibration. The actual numerical calculations for 
the transformed second mode as well as for the first mode and trans- 
formed third mode of an airplane wing vibrating at zero airspeed are 
discussed subsequently as a special case of flutter. 

The more general equations of airplane flutter at nonzero airspeed 
may be interpreted in such a way that they can be solved by a process 
analogous to that just described for coupled natural, vibration. 


APPLICATION OF THE ITERATIVE TRANSFORMATION 
METHOD TO FLUTTER 
Flutter Equations 


The differential equations of equilibrium for 
simple harmonic motion are 

a wing executing 

;H l ♦“>}§- v ♦ V 

(15) 

' S “ V + V* 

(16) 

These equations govern a motion represented by 


Y(x,t) = y(x)e imt 

(17) 

4 >(x,t) = 0 (x)e l£Ut 

(18) 


The use of the structural -damping coefficients gy and g 0 in equa- 
tions (15) and ( 16 ) is discussed in appendix B. The expressions 
PyY + P 00 and Q^y + Q 00 are the intensities of applied force and 

torque, respectively. For aerodynamic and inertia' forces and torques 
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due to air flow and distributed mass, the P and Q coefficients have 
values given by the following formulas (rearranged from those in refer- 
ence 4) : For P , 

y 



P 0 = P R0 iP I0 

in' which 



in which 


Q Ry " 



and 


^Iy 


/b\ 3 r (l \2f|/7\ 

y \y + (w 


r 7 \ , 2 

-) b.co 


(19) 


( 20 ) 


( 21 ) 


( 22 ) 


(23) 


(24) 


(25) 


(26) 


(27) 
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And for Q^, 


in which 


= %0 ' iQ I0 


A / b \T A 2 \ 2G A \ 2F 1 2 2 I fy \ . £ 

Q *0 = fc) Q* ' > T I 5 + a j^ + 8 + a + " r j(iTj 0 b o 


*>< • ft)*® ■■*•(*• -)f - (s - ■*)f](5). 




For inertia forces and torques due to concentrated mass, the intensities 
of force and torque are, respectively, 

, Ky + P00 , . 

v + v* - y 1 < 31 > 

J y dx— >0 


in which 


= lim 
dx— >0 


■ ii) b 

{if b &Y4)o b °“ 


y - (tf (“ k A 1 
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For a cantilever wing the boundary conditions on the displacements 

are 


II 

15 

11 

? 

II 

O 

II 

El (l + ig y ) 

/ x =0 

_ dx^ 

S EI ( 1 + 1^^1 T ■ 

1 dx_Jx=L L 

wf 1 * ig 0 )^" 

— Jx= 


(36) 


The differential equations (15) and (l 6 ) are now written with the 

eigenvalue <sr as an explicit factor. Thus equations ( 15 ) and (16) 
become 

2 2 

^2 E^l + igy )^ = o) 2 (P y 'y + p 0’0) (37). 

and 

■ s ^i 1 + ig 0)i! = ^ (V y + y^) (3 8 ) 

in which the P' and Q' coefficients are equal, respectively, to 
the P and Q coefficients divided by to 2 . 


Formulation of Pseudoflutter Problem 

Thoje solutions o> 2 , (y,0) of equations ( 37 ) and ( 38 ) for 

which a) is a real and positive (not complex) quantity represent the 
steady harmonic motions of true flutter. However, because the P and 
Q coefficients are in general complex and because of the presence of 
structural damping, the solutions of equations ( 37 ) and ( 38 ) will, in 
general, be complex and will include complex eigenvalues or. As in 
other methods of flutter analysis, the problem is made tractable by 

assuming at the beginning a value of the parameter k = This assump- 

. tion fixes the values of the P and Q coefficients. A real value 
of k is assumed because v must be real and only real values of cd 
can represent flutter. To avoid the inconsistence of assumed real 
values of k and obtained complex values of in the solutions, the 

problem itself is altered by introducing an artificial damping so that the 

m 2 

complex eigenvalue is given by ■ 1 ~ . , where g is the coefficient 

± + 1 S a a 

of artificial damping. Thus the differential equations of what may be 
termed the pseudoflutter problem become 
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<3x2 EI ( 1 + lg y)^2 1 + i8 n ( P y' y + P 0'0) 


dx c 


(39) 


S “l 1 + ig 0)to ‘ r^ir^y' y + V*) 


(40) 


p 

The value of a> can now be real for any assumed real value of k and 
is therefore the square of the frequency of the steady harmonic motion 
maintained by the artificial -damping forces and the naturally present 
aerodynamic, inertia, structural, and structural -damping forces. True 
flutter is possible for those special cases in which g a is zero. 

Equations (39) and (40) are similar in form to equations (7) and (8) 
and can be solved by the iterative transformation procedure in a way 
completely analogous to the solution of the ordinary problem. The com- 
plications introduced by the presence of air forces require, however, 
that a set of solutions be obtained for each of several assumed values 
of k. The fact that most of the /functions involved are complex virtu- 
ally quadruples the labor as compared with that required in the ordinary 
coupled natural -vibration problem. 


Steps in the Iteration as Applied to Flutter 

The iteration procedure employs the basic differential equations (39) 
and (40) in their integral forms which, for the cantilever wing under 
consideration, are 


y 




+ P0 , 0)(dx) 1 *' 


(41) 



1 



(Qy'y + Q0'0)(dx) 2 


(42) 


1 + ig a 

in which C stands for the more convenient form p — of the 

Qj 

eigenvalue. The iteration of equations (4l) and (42) follows the same 
form as the iteration of equations ( 9 ) and (10). Briefly, the steps 
are as follows: 
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(1) A real value of k is assumed and the values of the complex P 
and Q coefficients are computed. 

(2) An assumption is made for the desired mode y,0. (in the first 
cycle of iteration the assumed mode may be real but in the following 
cycles it will be complex.) 


(3) The complex loadings 



and Q^y + Q00 


are computed. 


(4) The integrations indicated in the equations are carried out 
numerically to get the complex intermediate derived mode. 


(5) The sweeping operations are performed by using the complex 
lover-order transformed modes previously determined. For convenience 
in numerical calculations, the flexural and torsional components of the 
complex derived (swept) transformed mode are computed in the forms ' 


1 7 ° L h 

C Vo C 


■*y 


(43) 


and 


1 lA ®o*o „ 

c Vo G o J o L 2 ^ 


(44) 


respectively, in which K and 
of the spanwise coordinate x. 


YLj are nondimensional complex functions 


(6) The derived and assumed modes are compared by computing their 
ratios at the stations of the wing. If these ratios are not reasonably 
the same, additional cycles of iteration are carried out until the ratios 
are reasonably the same. In the limit (never obtained in practice) the 
ratios will be identical and the proper value of C is that value which 
makes them unity; that is. 


1 7 o iA . K 1 7 o Vo b o 
C u Q E 0 I 0 oV = C n Q E 0 I 0 GqJ q l 2 *0 

y 0 


(45) 


in which y and 0 constitute the assumed mode of the limiting cycle 
and the functions in the numerators constitute the derived mode of the 
limiting cycle. 
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Equation (45) may be stated in the form 

in which D and H are nondimens ional real numbers. 
1 + ig a 

— , equation (46) may be written 


(46) 


Inasmuch as C 


is defined as 


CD 


(D + iH)— 


7q 


1 + igc 


^o ®o^o CD^ 

from which the frequency and artificial -damping coefficient are 
obtained as follows : 


(47) 


co = 




EqIo ^o 


L 4 J o 


(48) 


g 



(49) 


The relative air velocity corresponding to the assumed value of k is 
given through the definition of k, that is, 



(50) 


NUMERICAL EXAMPLES 


Numerical computations presented in this section illustrate the 
actual application of the iterative transformation procedure first to 
the ordinary natural -vibration problem (vibration at zero airspeed) and 
then to the flutter problem. All examples deal with the cantilever wing 
shown in figure 3- 

The geometric, structural, and mass properties of the wing are given 
in figure 3. A station coordinate system is employed for the purposes 
of the required numerical integrations. Four stations along the span 
have been selected as indicated in the figure j one of these stations is 
located at the spanwise position of the concentrated mass. The 
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distributions of forces and displacements over the span are considered 
to be adequately defined (through interpolation) by the forces and dis- 
placements at the four selected stations. The selection of a system of 
stations in any problem is important because it greatly influences the 
amount and accuracy of the work to follow. In problems , such as the 
present one, that involve concentrated masses, a station must be placed 
at each concentrated mass because displacements at the concentrated 
masses must be known. (More generally, a station must be placed at each 
discontinuity. Discontinuities may be present in the distribution of 
the structural stiffness and in the plan form as well as in distribution 
of mass. ) The other stations should be equally spaced between the dis- 
continuities, and for the system of parabolic interpolation used in the 
numerical integrations in this paper there must be a minimum of one sta- 
tion between each adjacent pair of discontinuities. The total number 
of stations should be the smallest possible that is consistent with the 
desired accuracy because the calculation effort increases rapidly with 
an increasing number of stations. In coupled systems, the number of 
degrees of freedom allowed is twice the number of stations selected; 
that is, the number of degrees of freedom in either the flexural or 
torsional component of displacement is equal to the number of stations 
employed. Experience has indicated that with parabolic approximations 
results accurate to at least two significant figures in the highest mode 
computed can be obtained by employing numbers of stations as follows: 

For uncoupled systems, the number of stations should be two greater than 
the order of the highest mode to be computed; for coupled systems, the 
number of stations should be one greater than the order of the highest 
mode to be computed, with a minimum of three stations. More than these 
minimum numbers of stations may be required if their use is dictated by 
sufficiently many discontinuities. 


Ordinary Coupled Natural Modes and Frequencies 

The calculations for the first, second, and third modes at zero 
airspeed for the wing of figure 3 are shown in tables 1, 2, and 3 , 
respectively. In this case k = 00 and the only aerodynamic forces are 
the apparent-mass forces. For simplicity, structural damping is dis- 
regarded; therefore, all quantities entering the problem are real. The 
numerical values of the aerodynamic-inertia force coefficients for k = <*>, 
as well as for other values of k to be used subsequently, are given in 
table k. 

The first coupled mode is computed in table 1. Table 1 shows in 
separate tabulations the flexural and torsional parts of the calculation. 
The first cycle of iteration (part (a) of the table) is shown in full 
detail. Two forms for the torsional part of the calculation are shown: 

The first form may be used when the torsional stiffness GJ is constant 
over each bay or over the whole length of the wing; the second form. 
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which requires slightly more work, must he used when GJ is variable 
and may be used, as in this case, when GJ is constant. The second 
and third cycles of iteration are summarized in parts (b) and (c) of 
table 1. 


Details of the first cycle of iteration, if the procedure that 
applies only for constant torsional stiffness GJ for the torsional 
part of the calculation is used, are as follows: In columns 1 of 


table l(a) the two parts 


(1) 

yq and 



of the assumed first mode are 


listed. The torsional component is assumed to be zero because it wili 
ultimately be small and is difficult to estimate. Columns 2 and 3 are 
the appropriate products of the assumed mode and the distributed-force 
coefficients. Columns 4, which are the sums of columns 2 and 3> give 
the two components of the external load which correspond to the assumed 
mode and the arbitrary frequency od. Columns 5 give the concentrated 
loads (external forces and torques) that are equivalent to the distrib- 
uted loads of columns 4. These equivalent concentrations are given in 
columns 5 in terms of the pertinent distances between stations Xq and 
in columns 6 in terms of the reference distance X Q . Formulas used for 
computing the equivalent concentrations from the distributed loads are 
given in appendix C. Columns 7 and 8 are the appropriate products of 
the assumed mode and the concentrated-force coefficients. Columns 9 
are the total concentrated loads, the sums of columns 6, 7, and 8. 


The flexural and torsional calculations must now be described 
separately. In column 10’ for flexure, the average shears in the bays 
between stations are found by a successive summation of the concentrated 
loads from the tip where the shear is zero inboard to the root. In 
column 11 the increments of bending moment are computed by multiplying 
the shears by the bay lengths in terms of X Q . The bending moments of 
column 12 are found by a successive summation of the increments of 
bending moment from the tip where the bending moment is zero inboard to 
the root. Column 13 gives the distribution of curvature, which is 
obtained by dividing each ordinate of the bending-moment curve by the 
local value of El (El in this example is constant). Equivalent con- 
centrated curvatures are now obtained by applying to the distributed 
curvatures the previously used formulas for equivalent concentrations. ' 
Column 14 gives these equivalent concentrations in terms of the dis- 
tances . Xq, and column 15 gives them in terms of the reference dis- 
tance X Q . The average slopes in the bays are obtained in column 16 by 
a successive summation of the concentrated curvatures from the root where 
the slope is zero outboard to the tip. The increments of derived 
flexural displacement are computed in column 17 by multiplying the aver- 
age slopes by the bay lengths in terms of X Q . The flexural compo- 
f 2 V 

nent y^ > of the derived mode is obtained in column 18 by a successive 
summation of the increments of displacement from the root where the 
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displacement is zero outboard to the tip. Column 19 gives the ratios 
at the selected stations of the derived flexural component to the 
assumed flexural component. 

Columns 10 to 15 for torsion are now considered. Column 10 gives 
the average twisting moments in the bays of the wing and is obtained by 
a successive summation of the concentrated torques of column 9 from the 
tip where the twisting moment is zero inboard to the root. The average 
twists in the bays are computed in column 11 by dividing the average 
twisting moment in each bay by the local value of GJ (GJ in this 
example is constant over the whole span). The increments of derived 
torsional displacement are obtained in column 12 by multiplying the 
average twists by the bay lengths in terms of X 0 . The torsional compo- 
nent of the derived mode is computed in column 13 by a successive sum- ' 
mation of the increments of displacement from the root where the dis- 
placement is zero outboard to the tip. Inasmuch as the derived displace- 
ment of column 13 is in terms of GJ, the displacement is converted 
into terms of El in column 14 so that it may be compared with the 
assumed torsional displacement on the same basis as the assumed and 
derived flexural displacements are compared and so that the next cycle 
may be started with displacement components having the same dimensions 
as the assumed mode of this first cycle. Column 15 normally would con- 
tain the ratios at the selected stations of the derived torsional compo- 
nent to the assumed torsional component, but in the case of table 1 (a) 
these ratios are meaningless because the torsional component will ulti- 
mately be different than was assumed in column 1 . 

Before the results of further cycles of iteration for the first 
mode are described, the form that the numerical integration for the 
torsional component must take when GJ is variable is described. In 
the part of table 1(a) showing the calculation for variable GJ, 
columns 1 to 4 are the same as in the calculation for constant GJ. The 
form of the numerical integration changes at column 5 . Column 5 consists 
of increments of twisting moment over the bays. These increments are 
obtained as increments of area beneath the curve of distributed torque 
(column 4). Formulas used for computing these increments are given in 
appendix C. In column 5 the increments of twisting moment are given in 
terms of the distances X^, and in column 6 they are given in terms of 
the reference distance X 0 . The twisting moments at the selected sta- 
tions due to the distributed torsional loading are obtained in column 7 
by a successive summation of the increments of twisting moment. The 
components of external concentrated torque are -as for constant GJ and 
are given in columns 8 and 9* The applied concentrated torque gives 
twisting moments as shown in column 10. Column 11 is the sum of 
columns 7 and 10 and gives the total twisting moments at the selected 
stations. (Note that in columns 10 and 11 there is a discontinuity in 
twisting moment at the station having the mass discontinuity.) Column 12 
gives the distribution of twist found by dividing column 11 by the local 
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value of GJ (GJ being in general not constant). The increments of 
derived torsional displacement are computed in columns 13 and 14 by 
applying to the values of column 12 the same formulas applied previously 
to column 4. The torsional component of the derived mode (columns 15 
and 16 ) is, except for small computational discrepancies, the same as 
in the previous method, as it should be. 


Two additional cycles of iteration were found to be adequate for 
the determination of the first mode and frequency. The results of 
these iterations are shown in parts (b) and (c) of table 1. In table 1(b), 
for example, columns 1 give the two components of the assumed mode of 
the second cycle, which are obtained by normalizing the derived mode of 
the first cycle to unity in the flexural component at the tip station. 

This normalization is not essential but facilitates manipulations and 
comparisons by keeping the numerical values in all cycles within the 
same range of magnitude. Columns 2 give the derived mode obtained by 
the numerical integration procedure just described. The ratios of 
derived to assumed mode are given in columns 3 for both components of 
displacement. These ratios are seen to be fairly uniform. The ratios 
obtained in the third cycle in table l(c) axe, for practical purposes, 
identical. The averaging device shown in columns 4 of table 1(c) and 
to the right of table 1(c) is adopted as a quick and generally quite 
accurate way of smoothing out small discrepancies that remain in the 
ratios after convergence is almost complete. This device, although 
clearly not necessary in the case of table 1(c), is useful in other 
cases throughout the numerical examples and is explained as follows: 

The two ratios in columns 4 are obtained by considering the flexural and 
torsional components of the displacement separately and then dividing 
the sum of the station values of the derived displacement by the sum of 
the station values of the assumed displacement. When a discrepancy 
remains between two ratios of the type in columns 4, the average of these 
two is taken as the final value; the final value for this case is given 
in the calculation to the right of table 1(c). This device gives greater 
weight to the larger ordinates and is in that respect similar to other 
weighting procedures such as the energy and least-squares methods but is 
much simpler. If the assumed and derived displacements contain both 
positive and negative ordinates, the negative ordinates should be changed 
to positive for the purpose of the summations. The remaining calculation 
shown to the right of table 1(c) gives that value of the arbitrary fre- 
quency which makes the ratio just computed unity. As proved in 
appendix A, this value of o> is the fundamental frequency U 3 q. 


Table 2 gives the main results of three cycles of iteration 
required to obtain satisfactory approximations of the second frequency 
and the transformed second mode at zero airspeed (k = »). Columns 1 of 
the first cycle (parts (a) of table 2) contain the two components y(p ' 

and 0^2^ assumed transformed second mode. This mode must have 
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one zero ordinate (excluding the root ordinates). Although this zero 
ordinate may theoretically be taken at any station, the numerical 
accuracy of the results is greatest if the zero ordinate is placed at 
the station and in the component where the preceding mode (the first) 
has its maximum numerical value (since the numerical process is such 
that the larger ordinates contain more significant figures than the 
smaller ordinates). Therefore, the zero ordinate of the transformed 
second mode is placed at the tip station in the flexural component, this 
location being designated station A. In the numerical values of 

columns 1, the flexural component y^^ w ould normally be taken as zero. 


(The values that are shown are estimated from flutter calculations that 
had previously been made for this wing.) Columns 2 give the intermediate 
derived mode obtained by numerical integration. Columns 3 constitute 
the first -mode sweeping function. The shape of this sweeping function 
is given by columns 2 of table l(c) and its magnitude is such as to be 
equal and opposite to the intermediate derived mode at station A. Thus 
the derived transformed second mode (columns 4), which is the sum of 
columns 2 and 3, has zero ordinate in the flexural component at sta- 
tion A and a shape comparable to the assumed mode, as indicated by the 
ratios -in columns 5. The ratios in the next two cycles (parts (b) and 
(c) of table 2) show marked improvement in uniformity. The final value 
of the ratio computed below the table gives, as proved in appendix A, the 
value of the second frequency uq , as shown. 


The main results of the iterations to obtain satisfactory approxi- 
mations of the third frequency and the transformed third mode at zero 
airspeed are stated in table 3. Typical operations required in a cycle 
are outlined in table 3(a). Columns 1 give the assumed transformed third 

mode made up of the two components y^^ 8X1(1 0^^ . transforme(1 

third mode is to have a zero ordinate in the flexural component at the 
tip station as in the transformed second mode and a zero ordinate in 
the torsional component at the tip station. The location of the second 
zero ordinate is designated station B. To -obtain greatest numerical 
accuracy, the selection of the second zero ordinate is governed by the 
same rule that was used for selecting the first zero ordinate, namely, 
that the new zero ordinate should be placed at the station and in the 
component where the preceding mode (the transformed second) has its 
maximum numerical value. The numerical values that are shown in 
columns 1 are estimated from previous flutter calculations: the torsional 

component 0^ would normally be taken as zero. Columns 2 give the 

intermediate derived mode, and columns 3 give the first-mode sweeping 
function which, as before, has a magnitude at station A that is equal 
and opposite to the intermediate derived mode. Columns 4 constitute the 
transformed-second-mode sweeping function which has a shape given by 
columns 4 of table 2(c) and a magnitude at station B equal and opposite 
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to the sum of the intermediate derived mode and the first -mode sweeping 
function (the sum of columns 2 and 3). The derived transformed third 
mode of the first cycle is the sum of columns 2, 3> and 4 and is given 
in columns 5* The ratios in columns 6 are far from uniform. The ratios 
in the second and third cycles (parts (b) and (c) of table 3) show 
improvements in uniformity. The iteration is discontinued at the end of 
the third cycle where the ratios are about as uniform as they can get 
with the limited number of significant figures that are present. The 
frequency obtained by the smoothing device is the third frequency 
and has the value shown. 

The patterns laid out in the foregoing examples establish the 
general technique that can be used to obtain zero-airspeed modes' and 
frequencies higher than the third. Guiding rules for determining the 
number of selected stations to be employed have been given previously. 
These examples also set the basic pattern for the computation of the 
modes and eigenvalues of pseudoflutter and of flutter. 


Modes and Eigenvalues of Pseudoflutter and of Flutter 

The operational solution in reference 5 gave for the wing under 
consideration (fig. 3) a reduced frequency at flutter of 0.1443. In 
order to use this operational solution, this same value (k = 0.1443) is 
used in the flutter calculations that follow. 

The calculations for the first, second, and third modes at 
k = 0.1443 are shown in tables 5* 6, and 7* respectively. Aerodynamic- 
inertia force coefficients have been computed by equations (19) to (35) 
and their values are given in table 4. Structural damping is disregarded, 
although a note on the method of incorporating structural damping in 
the calculations is made subsequently. 

Table 5(a) shows in detail the first cycle of iteration for the 
first mode. The form of the computations is the same as that shown 
previously for the determination of zero-airspeed modes. The amount of 
computation, however, is between three and four times that required for 
zero-airspeed modes because of the fact that the functions involved are 
complex and thus must be described by two parts - a real part and an 
imaginary part. Columns 1 and 2 are the real and imaginary parts, 
respectively, of the assumed first mode. As a start, all parts of the 
assumed mode except the real part of the flexural component are taken 
as zero.. Columns 3 to 6 are the real parts of the products of 
aerodynamic-inertia coefficients and the assumed mode, and thus their 
sums (columns 7) are the real parts of the distributed load. If the 
expressions for the distributed load are considered, this condition is 
more evident. The distributed forces producing flexure are given by 
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( P Ry - lP iy) ( y R + ly i) * ( P R0 - 1P I0) fa + **i) * Vr + p r/r + 

P Iy y i * P I0^I + 1 ( P Ry y i + P R0^I " P Iy y R ‘ P I0 ^r) 

The terms of the real part of equation (51) appear in columns 3 to 6 in 
the flexural part of table 5(a); the terms of the imaginary part of 
equation (51) appear in columns 22 to 25 in the flexural part of 
table 5(a). This separation of real and imaginary parts allows the dis- 
placement due to each part to be computed separately. A similar explana- 
tion can be made for the quantities in columns 3 to 6 and columns 18 
to 21 in the torsional part of table 5(a). 


Real and imaginary parts of the concentrated loads that are equiva- 
lent to the distributed loads are computed as explained previously by 
the formulas of appendix C. These values are shown in columns 8, 9, 27, 
and 28 in the flexural part and in columns 8, 9, 23, and 24 in the 
torsional part. The real and imaginary parts of the loads due to the 
concentrated mass follow next in order, and the total concentrated loads 
are given in columns 12 and 31 in the flexural part and in columns 12 
and 27 in the torsional part. The average shears, average twisting 
moments, and bending moments are then computed as described previously. 


The remaining parts of the computations in table 5(a) that ar ( e 
associated with the real parts of the load are described as follows 
(the remaining parts that are associated with the imaginary parts of the 
load are similar): Column l6 in the flexural part gives the distributed 

curvature due to the real part of the load. This curvature is obtained 
by dividing the ordinates of the real part of the bending -moment curve 
by the local values of the complex flexural stiffness EI^l + i %)( 1 + ig a 

In these examples, any actual structural damping is disregarded; there- 
fore gy is zero. The factor 1 + ig a , containing the as yet unknown 
artificial-damping coefficient, combines with o 2 to give the factor l/c 

in column l6, C being the arbitrary eigenvalue — If the actual 

ur 

structural damping gy is regarded as other than zero, the values in 
column 16 would be computed as follows: The real and imaginary parts of 

the bending moment would be combined into the complex bending moment 
Mr + iMp This complex bending -moment distribution would then be divided 


Mp I Xl'lT 

by the local values of the complex stiffness to give — -p: err— p. 

l 1 + ig y)( + ig a) 

The factor 1 + ig a would be carried along in the arbitrary eigenvalue C, 

Mr* + iMT 

and the numerical values of the real part of the quotient — $ 

E1 1 * lg y) 


+ liVh 
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would be placed in column 16. The imaginary part of the quotient would 
be similarly placed (in column 35) in the calculations associated with 
the imaginary part of the load. The average twists due to the real part 
of the load are computed in column 14 in the torsional part of table 5(a), 
and those due to the imaginary part of the load must also be computed. 
These calculations follow the same pattern as those just explained for 
the curvatures. The complex torsional stiffness Gj(l + ig^) (l + ig a ) 

enters in place of the complex flexural stiffness. If GJ or gj is 
variable over a bay length or over the whole span, the numerical ^ 
integration for the torsional part of the calculations should be carried 
out as explained in the part of table 1(a) that deals with variable GJ. 


The numerical integrations are completed in the manner already 
described, and the derived mode is thereby obtained in the form of four 

components of displacement. The flexural components are y^ and y^ 

t UR IX 

of columns 21 and 40 in the flexural part. The torsional components 

axe 0^ and 0^ of columns 17 and 32 in the torsional part. How- 
ever, these components are not actually the real and imaginary parts of 
the flexural and torsional components of the derived mode, because each 
one of them contains the complex factor 1 + ig a . Nevertheless, the 

complex derived mode is given by yjj^ + iy^ and 0^ + i0^ . 

The complex ratios of the complex derived mode to the complex 
assumed mode are computed in column 4l in the flexural part and 
column 33 in the torsional part. Only two of these ratios have actually 
been computed but they are sufficient to indicate the need for further 
cycles of iteration. 


A total of four cycles of iteration (the main results of the last 
three are shown in parts (b), (c), and (d) of table 5) were required for 
satisfactory convergence. In columns 6 of table 5(d) and immediately 
below table 5(d), the smoothing device described previously is applied 
to obtain the best single value of the ratios. The fundamental (first) 
eigenvalue is that val^e of C which makes the ratio unity. Thus 

Ci = (269.5 - 82. 2i] , and since C 1 is defined as the 

t ° 1 

frequency and artificial damping of the first mode are obtained from 
the real and imaginary parts of the equation 


1 + ig , ! A 

■ ( 269 - 5 - 82 - 21 )*t 


(52) 
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The calculation of these quantities and the corresponding airspeed v^ 

htni 

which is obtained from the relation v^_ = are shown at the bottom 

of table 5* 


Tables 6 and 7 show the main results of the iterations to obtain 
the transformed second and third modes for k = 0.1443. Four cycles of 
iteration for each mode gave satisfactory convergence. The assumed 
modes of columns 1 and 2 of tables 6(a) and 7(a) were taken in the forms 
recommended previously in connection with tables 2 and 3* In tables 6 


(n) 


and 7, the complex intermediate derived modes are given by y^ R 
<(a) ^ i0(n) 


+ iy 


and 0. 

y (n) 

y blR 


bR 

+ iy 


(n) 

bl 


(n) 

bll 


bl ' 
and 


the complex first-mode sweeping functions, by 


0 - 


(n) 

blR 


10 


(n) 

bll 


with shapes corresponding to columns 3 


and 4 of table 5(d); and the complex transformed-second-mode sweeping 
functions, by y^ p + iy^ T and 0^ R +■ 10^1 with sha P es corre- 


'ba2R ,7 ba2I “ ^ba2R 
sponding to columns 7 and 8 of table 6(d) 

below table 7(d) give for the third eigenvalue, g = 0.030 

a J 

03 = 168.9 radians per second. The corresponding airspeed is 
v 3 = 390 feet per second. 


The results computed in and 
and 


Computation of True Modes 

Because the critical flutter velocities are given directly by the 
eigenvalues, knowledge of the true modes in flutter problems is of no 
value (at least of no value recognized at present). The same statement 
applies to the transformed pseudoflutter modes, with the exception that 
in the iterative method their determination is a necessary adjunct to 
the determination of the eigenvalues. In ordinary problems of forced 
vibration (at zero airspeed), however, the true modes are often used 
with great advantage. For this reason and for the sake of completeness 
of the presentation of the iterative transformation procedure, the 
method of determining true modes from results of the iterative trans- 
formation procedure is illustrated in tables 8 and 9- 

The computations in tables 8 and 9 pertain to the same wing analyzed 
in the previous examples. The modes computed are for k = 0.1443. The 
true third mode as computed in table 9 may therefore be compared with the 
flutter mode computed for this wing by the operational method in 
reference 5. „ 

In table 8, the true second mode is computed as follows from func- 
tions appearing in the last cycle of iteration for the transformed 
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second mode (table 6(d)): Preceding the table proper is the calculation 

C 1 

of the eigenvalue factor F-^ = — - 1 that is needed for computing the 

true second mode. In the terminology of tables 6(d) and 8 and as shown 
in appendix A, the true -second -mode shape is given by 


and 


in which 


and 



+ 


10 


< (4) 

dIR 


+ l' 


i0. 


W 

bll 


II 


Cl- 

Co 


- 1 


(53) 


(5*0 


(55) 


(56) 


Columns 1 of table 8 show the key ordinate |y),i + iy!,^ of the 

rblR •'bll/ A 

first-mode sweeping function y^ + iy^h 0^ + i0^) as given in 

D1K Dll Dli\ bll 

columns 5 and 6 of table 6(d). The key ordinate is taken as the largest 
ordinate (the ordinate at station A) for the reason of accuracy 
cited previously. The key ordinate of the first-mode shape 

y + iy ,0, +10 (equal to the first terms on the right-hand sides 

1R II 1R II 

of equations (53) and (5*0) is shown as the boxed value in column 2 and is 
obtained by dividing the value in column 1 by the eigenvalue factor F 12 . 

The other values in columns 2 are obtained by using the key ordinate in 
conjunction with the first-mode shape given in columns 3 'and 4 of 
table 5(d). Columns 3 show the transformed-second-mode shape 

r (- \ / R \ / J- \ / r \ 

y a2R + ^a2R + ^a2I ( e< l ua i ihe second terms on the right-hand 
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sides of equations (53) and (54)) given by columns 7 and 8 of table 6(d). 
The sum of columns 2 and 3 which is given in columns 4 gives the shape 
of the true second mode y^ + iy^, ^2R + i ^2I ( e< l ua d to the left-hand 

sides of equations (53) and (54)). 


In table 9 } "the computation of the true third mode proceeds as 


follows: The necessary eigenvalue factors F-^ 


= C3 " 1 F 23 = cf * 1 

are computed as shown. In the terminology of tables 7(d) and 9 and as 
shown in appendix A, the true -third-mode shape is given by 


3R 



( y llR + iy linj + i y 12R 



+ 



+ iy 


(5) 

a3I 


(57) 


and 


^3H + i ^3I “Kir + ^lll) + (^12R + i « ZS 12l) + (*a2R + ^a2l) + (^3R + ^3I y 

(58) 

in which 

(f| " ^ ( y llK + Iy lll) * (gf * ^(^lSH * ly 12l) = y bIR + 1 ^‘ ,) 


'bll 


(59) 


(c| ' X ) (^11R + ^lll) + (c| ' ■*■) (^12R + i ^12l) = s 4 lR + i 4 ll 


(60) 


y a2R + iy a2I 


y ba2R + iy ^I 

c 3 


(61) 


0 Q on + 


a2R ^a2I 


0< 4 ) + i^ 4 ) 

^ba2R 1 ’°ba2I 


^ - 1 


(62) 
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811,1 y 12R + iy 12I' ^12R + ^121 1S t0 y a2R + iy a2I> ^a2R + i ^a2I as 

y lR + iy lF ^1R + ^11 1S to y a2R + iy a2I' ^a2R + 1 ^a2I in table ®* 

The key ordinates (y<^ * ijr<J>) A and (0^ + i0^i) B of the first 

and second sweeping functions appear in col umn s 1 and 2 and are taken 
from columns 5> 6, 1 > and 8 of table 7(d). The key ordinate of the 
functions ^y 12R + iy 12I ) , F^ (0 12R + i0 12I ) , which are equal to the 

second terms on the left-hand sides of equations (59) and (60 ), is com- 
puted in columns 3 by using the key ordinate of columns 2 in conjunction 
with ordinates at stations A and B in columns 2 and 3 of table 8 as 
follows : 




9 


( y LR + iy ll) A 

p) ' p)\ 

^a2R + 1 ^a2l/ 


JLltable 8 




table 


9 


(63) 


The key ordinate of the functions F^y^ + iy i;LI ), f i3(0hk + i0 11]; ) > 

which are equal to the first terms on the left-hand sides of equations (59) 
and ( 60 ), is given in columns 4 and, in accordance with equations (59) 

and ( 60 ), is the difference between y^l^ + iy^, 0^^ + i 0tili of 

columns 1 and F 2 3( y 12R + iy i2l ) * F 23(^12R + i ^12l) of columns 3- The 
key ordinates of the first-mode shapes y 11R + iy i;LI > 0 11R + i0 113 . and 

yi 2 R + iyi2I> 012R + -^0121 0X6 s ^ own in columns 6 and 5 and axe obtained 

by dividing the values in columns 3 and 4 by the appropriate eigenvalue 
factors. The sum of the key ordinates of columns 5 and 6, shown as the 
boxed value in column 7* is the key ordinate of the total-first-mode 
shape y-^ + iy^_j, 0 ir + i0u which is equal to the sums of the first ■ 

two terms on the right-hand sides of equations (57) and (58). The other 
values in columns 7 are obtained by using the key ordinate in conjunction 
with the first-mode shape given in columns 3 and 4 of table 5(d). The 
key ordinate of the transformed- second -mode shape y^ pR + iy^j, 0a2R + ^a2V 

which is equal to the third terms on the right-hand sides of equa- 
tions (57) and ( 58 ), is shown as the boxed value in columns 8 
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is obtained by dividing the value in colujnns 2 by the eigenvalue fac- 
tor F 23 . The other values in columns 8 are computed by using the key 

ordinate in conjunction with the transformed-second-mode shape given in 
columns 7 and 8 of table 6(d). Columns 9 show the transformed-third- 

mode shape y a ^p + iy a3 j , 0 a 3R + 10a3I ( equal to the fourth terms on 

the right-hand sides of equations (57) and ( 58 )) given by columns 9 
and 10 of table 7(d). The sum of columns 7 ) 8, and 9 given in columns 10 
gives the shape of the true third mode y^ + iy^j, 0^ + 10^1 (equal 

to the left-hand sides of equations (57) and (58)). 


DISCUSSION OF RESULTS 

Trends and Comparisons of Numerical Results 


Results of the computations shown in the preceding section of the 
paper together with results of similar computations based on other 
assumed values of k are given in figures 4 to 6. Figures 4 and 5 deal 
with the wing to which the concentrated mass is attached. Figure 6 
gives data of a similar nature for the same wing without the concentrated 
mass. The computed results obtained by the Rayleigh-Ritz and operational 
methods and the experimental results, all of which are given for this 
wing in references 5 and 6, are also recorded in figures 4 to 6. 

In part (a) of figure 4 the solid curves show the variation of the 
artificial -damping coefficient g a with airspeed in each of the first 
three solutions. For each assumed value of k a dashed curve is drawn 
through points that represent solutions for that value of k. Part (b) 
of figure 4 shows in a similar way the variation of the frequency cd 
with airspeed and the lines of constant values of k. The facts of 
particular interest that are shown by these plots are as follows: 

(1) The true flutter condition is given by the third solution for 
a value of k between 0.1443 and 0.1590 at an airspeed almost equal to 
that found in the experiment. Here the computed value of g a is zero. 

The computed frequency at true flutter is also in very close agreement 
with the experimental value. 

(2) The operational solution is in good agreement with the experi- 
mental solution, but the solutions obtained by the Rayleigh-Ritz method 
with three and four modes vary by 72 percent and 22 percent, respectively, 
from the solution obtained by the operational method. The operational 
solution is theoretically the most exact even though it involves summa- 
tions of finite numbers of terms of infinite series. However, as pointed 
out in reference 5* its use is limited in practice to wings of uniform 
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section. In the present example the results obtained by the iterative 
method would be expected to be better than the results obtained by the 
Rayleigh-Ritz method because the eight degrees of freedom used in the 
iterative method are much less restrictive than the three or four used 
in the Rayleigh-Ritz method. Although exact agreement of the results of 
any of the computational methods with the experimental results is not to 
be expected, the better agreement of the iterative solution as compared 
with the operational solution is at first surprising. On further 
observation, however, this agreement must be credited to a fortunate 
disposition of the errors involved in the iterative method because, in 
the case of figure 6, the relative order of agreement of the operational 
and iterative results with the experimental result is opposite to that 
in figure 4. 

(3) The trends of the solid curves representing the first and 
second solutions in figure 4(a) indicate that both may cross the zero 
artificial-damping axis at very large airspeeds. But this conjecture 
is of no practical interest so long as a curve (the third solution) that 
crosses at a lower airspeed exists. However, the question of whether 
the curve for some solution higher than the third could cross the zero 
artificial-damping axis at an airspeed lower than that at which the 
third solution crosses demands an answer. 

w Reasonable assurance that, among all possible solutions, the 
curve of third solutions in figure 4(a) crosses the zero artificial- 
damping axis at a lower airspeed .than any other is provided by the 
trends of the curves for constant values of k in parts (a) and,(b) of 
figure 4. The curves of k show that the curves representing the fourth 
solution will most assuredly lie above and to the right of the solid 
curves in figure 4(b) and probably below and to the right of the solid 
curve for the third solution in figure 4(a). The curves of k in fig- 
ure 4(b) are straight lines by definition ^k = Prediction of the 

courses of the curves of k in figure 4(a) cannot be made with much 
certainty. They have a strong tendency to proceed to the right, but it 
is easy to believe that upward or downward changes in their directions 
could take place. The curve for the fourth solution, however, would 
probably cross the zero damping axis at a value of v between 500 and 
600 feet per second in figure 4(a). 

Figure 5 shows and compares the amplitude and phase distributions of 
modes computed by the iterative transformation procedure and by the 
operational method for the wing with a concentrated mass. The first and 
second modes as well as the more important third mode from the iterative 
solution for k = 0.1443 are plotted, and the third mode from the 
iterative solution for k = 0.1590 is also plotted. The third modes 
from the iterative solutions for the two values of k . agree very well 
in shape with the flutter mode obtained in reference 5 by the operational 
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method, and the operational mode lies between the two iterative modes. 
Thus the agreement of the iterative and the operational methods is 
again evidenced. 

Figure 6 is a plot similar to figure 4 but relates to the behavior 
of the wing analyzed in figure 4 if the concentrated mass is not present. 
There is very little similarity in the data of the two figures. The 
most notable difference is that in figure 6 the true flutter mode appears 
in the second solution instead of the third as before and that the 
flutter speed is lower than before. Of interest is the occurrence of 
almost equal eigenvalues in the second and third solutions for k = 0.50. 
.The flutter speeds given in figure 6 by all methods of solution, 
including the Rayleigh-Ritz method, are seen to be in substantial 
agreement . 


CONCLUDING REMARKS 


The paper has described the iterative transformation method sug- 
gested by H. Wielandt and has demonstrated the use of the method in an 
orderly computation of critical flutter speeds. Numerical comparisons 
with solutions obtained by other methods and with experimental values 
have been made. The applications made in this paper show promise for 
future practical use of the method. 


Langley Aeronautical Laboratory 

National Advisory Committee for Aeronautics 
Langley Field, Va., January 17, 1951 
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APPENDIX A 

ON THE CONVERGENCE OF THE ITERATIVE TRANSFORMATION PROCEDURE 

Introduction 


The extensive existing literature on the eigenvalue problems is 
concerned almost exclusively with the class known as self-adjoint prob- 
lems^ in which the eigenfunctions and eigenvalues are real. In recent 
years, non-self -ad joint eigenvalue problems have received increasing 
attention. This class includes the flutter problem in which the eigen- 
functions and eigenvalues are generally complex. The literature 
referred to by Wielandt in reference 3 reveals that the non- self -adjoint 
eigenvalue problem and the transformation method for its solution have 
been given some attention since at least 1928. Wielandt' s own work 
constitutes probably the most extensive contribution on the subject. 

The discussion on convergence given herein is not contained in 
Wielandt ' s work and may be considered a rigorous proof if the following 
assumption is valid: that the equations (equations (4l) and (42)) for 

the system (the wing) under consideration have an infinite n umb er of 
solutions that form a complete set for any value of the reduced fre- 
quency k. In the subsequent demonstrations, the validity of expanding 
arbitrary displacement functions in infinite series of eigenfunctions 
depends upon the validity of the assumption. That complete sets of 
eigenfunctions do exist seems plausible enough to justify reliance in 
the conclusions. 


Basic Relations 

For any one of the true solutions of the eigenvalue problem, for 
example, the eigenvalue and eigenfunction equations (4l) 

and (42) may be written as 
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To make the notation more concise, let the coupled mode y m ,0 ffl be 
represented by v m . Then if y m ,0 m is substituted into the right-hand 

sides of equations (Al) and (A 2), the left-hand sides may be represented 
by C m w m . Furthermore, because of the linear character of the equations 
of the problem, substitution of the function series 


H a i w i 

i=l 


(A3) 


into the right-hand sides of equations (Al) and (A2) gives for the left- 
hand sides the function series 


00 



C i a i w i 


(A4) 


The coefficients a-^ are , in general, complex. The complex eigen- 
values C± are assumed in the subsequent proofs, except where stated 
otherwise, to be different from each other, and the eigenvalue having 
the largest modulus is defined as C^, the second largest, as C 2 , and 
so forth, so that 


> C, 


(A5) 


Expressions (A3) and (Ah) are the expansions, in terms of the 
eigenfunctions and the eigenvalues, of the functions previously referred 
to as the assumed and intermediate derived modes, respectively. The 
subsequent proofs of convergence are based upon the fundamental relation- 
ship that exists between expressions (A3) and (A4). 


Fundamental Mode 

The fundamental mode and eigenvalue are found by iteration 
according to the original Stodola procedure. In the present terminology 
and- notation, this procedure and its proof of convergence are as 
follows: The coupled mode assumed at the beginning of the first cycle 

of iteration in general contains some component of each of the eigen- 
functions; therefore its most general expression is 



& 1 " 1 


(A6) 
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The intermediate derived mode (which in this case is also the final 
derived mode inasmuch as no sweeping operation is required to obtain 
t the first mode) is for this first cycle of iteration 


( 1 ) ( 2 ) 

w b = W 1 = L_ c i a i w i (A7) 


The second and following cycles are begun with the final derived mode 
of each preceding cycle, and thus the assumed and derived modes of the 
nth cycle are 



oo 

r 

i=l 


_ n-1 
C i a i w i 


(A8) 



C 


n 

i a i W i 


(A9) 


In accordance with the definitions given in equation (A5), all 
terms on the right-hand sides of equations (a 8) and (A9) except the 
first are negligibly small in comparison with the first for large values 
of n. In the limit the fundamental mode is obtained as 


lim w[ n+1) = lim C 1 n a 1 w 1 ' (A10) 

n — ^ oo n — ^ oo 


and the fundamental eigenvalue is obtained from 


lim 
n — ^oo 


(n+l) 

v. 




, (All) 


Transformed Second Mode 

The initial assumption of the transformed second mode in general 
is of the form ^ * 


( 1 ) 

w a2 



b i 



(A12) 
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in which the arbitrary coefficients bi are in general complex and 
the subscript A refers to values of either the flexural or torsional 
components of the eigenfunctions at station A. More specifically, if, 
for example, the nodal (zero) point of w a2 se l e cted to be at sta- 
tion A in the flexural component, then the subscript A refers only to 
the flexural components of w^_, W 2 , w^, • . . and not to their tor- 

sional components. Thus each term of the series in equation (A12) sat- 
isfies the requirement that either the flexural or torsional component 
of the assumed mode be zero at station A. 

To simplify the subsequent work as much as possible, the eigen- 
functions are henceforth assumed to be normalized to unity at station A) 
thus 

(vi) A = 1 (i = 1, 2 , 3 , • • •) (A13) 

Equation (A12) now takes the simpler form 

- Wi) (A14) 

i=2 


The assumed mode given by equation (Al4) leads, according to equa- 
tions (A3) and (A4)) to the following intermediate derived mode: 

00 

w b 1) =Y1 b i( C i w i - C l v l) (A15) 

i=2 ' ' 

Sweeping of this intermediate derived mode with the first-mode shape 
(previously determined) leads to the derived transformed second mode of 
the first cycle as follows: 



(A16) 


When each succeeding cycle is begun with the derived transformed 
second mode of its preceding cycle , the various functions for the nth 
cycle are 


n) 

w a2 



(A17) 
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(n) n-1 / . 

w b = 1_ c i ( c i w i " c l w l) (A18) 

(n+1) V s - n , . 

w a2 = 2_ C i 'b i (w i - wj (A19) 

i=2 ' * 


The limits as n approaches infinity are 


and 


lim v^g +1 ^ = lim ( w 2 " w l) 

n — ^oo n — -^oo v 


(A20) 


4 


lim 
n — ^oo 


(n+l) 

w a2 


w 


a2 



( A21) 


Equations (A20) and (A21) show that convergence to the exact -transformed- 
second-mode shape w 2 - and to the exact second eigenvalue C 2 can 
he obtained theoretically. 


True Second Mode 


The key to computation of the true second mode is readily found in 

the simple case illustrated in figure 1. In this case the sweeping 

function of the final cycle of iteration would be the displacement 

/ 2 2 \ 

produced by the forcing load 7^og - “>l )J±, in which y 1 is the 

first -mode component of the transformed second mode y The sweeping 

function is designated by y^ which has a well-defined numerical 

value in the iteration. Thus the value of y-^ could be found from the 
equation 





(A22) 
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yi 


'bl 


(A23) 


a>i 2 


- 1 


The sum of y &2 given in the iteration and y n given by equation (A23) 


gives yg, the true second mode) that is. 


a2 


+ y x = y 2 - = y s 


(A24) 


By analogy, the true -second -mode shape in the general (complex) 
problem under consideration is found as follows: The limiting value of 

the sweeping function is, from equation (Al8), 


* w ( n )\ 

lim w- = lim ^ ' V b 2 U l 


>> . 

n — >oo n — >oo 


n/ c l 


(A25) 


The expression analogous to equation (A23) is 


n — >«> 


bl 


n 


Ei .1 


— lim 


n- 


(A 26) 


The expression analogous to equation (A24) is 


lim w, 
n — >oo 


, • (n) 

lim w ' ' 

(n+1) n — ^oo bl 


a2 


Ei . i 

C 2 


n 

= lim C 2 ^ 2^2 
n — >°° 


(A27). 


which gives the exact shape of the true second mode. 


Transformed Third Mode 

The first cycle of iteration for the transformed third mode begins 
with an assumed mode that has two zero values, one of these being in the 
same (flexural or torsional) component and at the same station (sta- 
tion A) as previously employed for the transformed second mode. The 
other zero value may be taken in the same component as was the first 
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zero value and at a different station (station B), or it may be taken 
in the other component at any station, including station A. Either of 
these possible selections for the location of the second zero value is 
indicated, in the following equations by use of the subscript B. The 
initially assumed transformed third mode may be written as 


,(U - 


y d. w, - v - — \ fw - w \ 

fcl 1 \ W 2 " W 1 /bV 2 


(A28) 


in which the arbitrary coefficients are complex. Each term of the 

series in equation (A28) is zero at station A by reason of the normaliza- 
tions stated in equation (A13), and each term is also zero at station B. 

The various displacement functions for the general (nth) cycle of 
iteration may be expressed as follows: The assumed mode is 


w£ J = Y1 C w. - w - (— (w - w ) 

a 3 fr 3 X; .1 \2 ~ W l/ B ' ^ YJ 


(A29) 


The intermediate derived mode is 


(n) 

w b 




c l w i 



(A 3 0) 


The result after sweeping the intermediate derived mode with a first - 
mode shape such as to make the sum zero at station A is as follows: 



Sweeping of the mode given by equation ( A 3 l) with a transformed- second- 
mode shape such as to make the sum zero in the flexural or torsional 
component (as the case may be) at station B gives the derived transformed 
third mode as follows: 
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(n+1) (n) 

w a 3 = w b ' 


4 n) 

Wi 


Wn 



/M\ 


(n) 

K 

1 w x 
A 

w b 

\"i ! 

w 2 

~ W 1 


— 


— 


( w 2 - w l) 


B 


i=3 


(^) B h - 


Ihe limits as n approaches infinity are 
(n+1) 


lim v, 
n — >00 


a3 


n 

= lim 
n — >00 


”3 - W 1 - (»2 - V l) 


and 


(n+1) 

v oO 

lim — /-V ■ = C. 
(n) 3 


n' — ?oo w 


a3 


(A 32 ) 


(A 33 ) 


(A 34 ) 


As shown by equations (A 33 ) and (A 34 ), convergence to the exact- 
transformed-third-mode shape w^ - w i - - vqj and to the 

exact third eigenvalue C3 can be obtained theoretically. 


True Third Mode 

Computation of the true third mode is explained by referring again 
to the simple problem of pure flexural vibration in which air forces 
are excluded. The transformed third mode in this simple problem would 
be given by 



y a 3 = y 3 ‘ y l 


(A35) 
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The total load required to hold the beam in equilibrium in the 
shape y &3 is ' 


2 2 
70> 3 y 3 - 70^ 


1 - 


( 7 1 - y± ) 
kT 0 - y n 


yi - 7<*2 


2/y3 - yi 


^2 


B 


(A36) 


If the beam is vibrating with shape y^ at frequency 003 , the inertia 
load is given by 


7a> 3 





(A37) 


The forcing load required is the difference between the total load 
(expression (A 36 )) and the inertia load (equation (A 37 )), that is. 



The displacement produced by this forcing load is 



(A39) 


and this displacement must be equal to the sum of the sweeping functions 
in the last cycle of iteration (if the iteration has been carried to com- 
plete convergence). The first sweeping function is of the first-mode 
shape and the second sweeping function is of the transformed-second-mode 
shape. If the expression (A39) is written in the form 
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each of the sweeping functions contained in the displacement produced 
by the forcing load is obvious. Thus 



in which y^ and y^p designate the first and second sweeping func- 
tions, respectively. Both of these functions have well-defined 
numerical values in the iteration. 


If now a simpler notation is adopted, equations (A4l) and (a 42) 
can be written as 


and 


in which 



(A43) 


(a44) 


(A45) 


'12 


_ ( i 3 ~ y i ^ 

^ - y i; 


B 


(a46) 
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and 


a2 



(A47) 


The true third mode is clearly given by the sum of equations (A35), 
(A45), (A46), and (A47 )j thus 

y 3 = y a3 + y ll * y 12 + y a2 < AWJ > 


The transformed third mode y^ is given directly in the iteration. 

The procedure for finding the other components on the right-hand side of 
equation (a 48) is as follows: Component y , by equation ( A44 ) , is 


y 

ba2 



(a49) 


Component y 12 is known when y^ is known because its relation to y a2 

was established previously in connection with the transformed-second- 
mode calculations (see equation (A24)). Component y is then found 

by equation (A43) as 


r U)~ 


'bl 


11 


CDo 2 

^2 - 1 


2 " 1 ) y 12 


(A50) 


By analogy with the foregoing case, the true third mode in the 
complex -eigenvalue problem is found as follows: The limiting value of 

the second sweeping function is (see equations (A31) and (A32)) 
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lim w> 


- lim 
n — > oo 


■ C >'R • h ■ <*» 


The limiting value of the first sweeping function is (see equations (A30) 
and (A31)) 


lim w. 


* lim yr w i 

n — \ 1 /A 


n — >oo 


lim C.V p - 1 - & - - 1 

l — >» 3 3 c 3 \ c 3 c 3/\ w 2 - w : 


3/\ 2 


(A52) 


The quantities analogous to y^_2 y»p of equations (AV 3 ) and (a44) 

are, for the present case, lim and lim v^. The latter 

n — > oo n — >oo a 

quantity is obtained from the relation analogous to equation (A^9) as 
follows: 


lim w 


lim w 


l 3 


= n— ^ <A53) 


The relationship of lim and lim is obtained from 

n — >oo ^ n — >oo 

equations (A26) and (A20) of the section dealing with the transformed 
second mode. Thus, 




(A5M 


The quantity analogous to y, , of equation (A4 3 ) is, for the present 
, . (n) 

case, lim w-^ and. is obtained by an equation analogous to equa- 
n — ^oo 

tion (A50) as follows: 


lim 
n — >00 


n 
W 11 


lim 
n ■ ^ 00 



lj lim w^ 
/ n — >00 


1 


= lim 

n — >00 


c 3 n d 3 




(A55) 


The exact shape of the true third mode w 3 is given by the sum of 
equations (A33), (A53), (A54), and (A55), which is 


lim 


n- 


( W < 


(n+l) 

»3 


(n) 

w a2 


(n) 

W 11 


(n)\ 

+ w 12 ) 


lim d^w^ 
n — >00 


(A56) 


Fourth and Higher Modes 

Extensions of the proofs to modes higher than the third can be made 
in a manner similar to the foregoing proofs. By this means, the itera- 
tive transformation procedure can be proved, under the assumptions 
stated at the beginning of this appendix, to be convergent for all modes 
and eigenvalues. 
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Cases of Eigenvalues Having Equal or 
Nearly Equal Moduli 

For a representative case, suppose that 


and that 


' or that 



(A57) 


(A58) 


(A59) 


Under conditions (A57) and either (A 58 ) or (A59), the assumed and 
derived modes after a few cycles of iteration will he virtually as 
follows (see equations (A17) and (A 19 )): 


“la* - c 2 n 'S(“2 - “]) + c 3 n ' 1 f 3 (w 3 - wj (A60) 

w a 2 +1 ^ ' C 2%(' , 2 - w l) + C 3 nb 3 ( u 3 " *]) < A6l > 

If |c 2 | is only slightly greater than (c^j, the second terms on the 

right-hand sides of equations (a 6 o) and (A 6 l) become negligibly small 
very slowly as n increases, even though they do become negligibly small 
as n approaches infinity. If |c 2 | and ) C 3 1 are equal, these terms 

never become negligibly small. .Thus, the problem of circumventing this 
slow convergence or apparent lack of convergence arises. 

A satisfactory method for coping with these conditions is to 
combine linearly the results of the last two cycles of the series of 
iteration cycles that have been performed. For best results in an 
actual problem, not less than the third and fourth cycles should be used 
for this purpose in order to reduce as much as practicable the effects 
of all higher -order components. 
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The following formulas for combining the results of the last two 
cycles are based on the assumption that the assumed and derived modes 
in each of the cycles contain only components of the types in equa- 
tions (A 60 ) and (A6l). 


The two components (with shapes w 2 - w^_ and w^ - w 1 ) clearly 

appear in the last cycle in proportions different than in the preceding 
cycle. (The proportion in each cycle is a complex function of the 
spanwise coordinate. ) Because of this differing proportionality the 
results of cycles n - 1 and n can be linearly combined so that the 
combined functions contain only one of the components W2 - wp 

and- W3 - wp. Accordingly , the ratios of both the flexural and tor- 
sional components of the combined functions at all stations should be 
equal to each other. In algebraic terms , this statement means t ha t 


' rw. 


1 rw 


(n) (n+lj\ 

a2 + w a2 
a2 + W a2 


= R 


(A62) 


in which r and R are (complex) constants, and the subscript S 
designates that the ratio may be evaluated at any station S , that is, 
that R has the same value for all stations. All w functions must 
be the same type of component, either flexural or torsional. 


Since S can be any station, the equality 


rw, 


(n) ( n+1 )' 


a2 


+ w, 


a2 


irw 


(n- 1 ) ^ J n ) 


a2 


+ w 


a2 


L 


(n) 

rw a2 


(n+1)' 
+ w a2 




,(n-l) ^ „(n) 


a2 


+ w' 


a2 


(A 63 ) 


exists, in which stations 1 and 2 must be different or may be the same, 
depending on whether the w functions on the left-hand side are the 
same or different types of components than those on the right-hand side. 
The two values of r that satisfy equation (A63) are 


A (n-l),(n+l) //a( n_1 > > ( n 1 A 2 A (n), (n+1) 

T^TTW 4 JhTn-l),(n) ] - ^:T),(n) < a61 *> 



NACA TN 2346 


49 


in which 


(n-l),(n) 


(n),(n+l) 




(n) 

r 

a2 


(n) 

/ 

a2 


w 


(•n 

r\ 


a2 /2 

w (n+1) ) 

. a2 j. 


(A 65 ) 


(a66) 


A (n-l),(n+l) = 


/ (n-l)\ , 

( n-l)\ 


(v ), 1 

\a2 j 

2 

L (n+1 ») | 

f (n+1)' 

W , 

) 

l® 2 k 1 

\ a2 t 

2 


The corresponding values of R axe 

(n-1), (n+1) // (n-l),(n+l)\" (n),(n+l) 


R = 


A (n-l),(n) yl A (n-l),(n) 


L (n-l),(n) 


(A67) 


(n),(n+l)/(n-l),(n) 

= (A68) 

These values of R are equal to and C ^ > and the corresponding 

values of r ; when placed in the expression 


rw( n ) + v( n+ -*-) 
a2 a 2 


(A 69 ) 
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give modes of the shapes wg - w-j_ and w^ - w^. When | | is nearly 

equal to | C 3 1 , the appropriate set of R and r to give the lower 
transformed mode Wg - is evident. When | Cg | and | | are equal, 

the mode obtained by equation (A 69 ) with either value of r may be used 
as the transformed second mode, but the trends of the eigenvalues that 
have been or will be determined at other values of the reduced fre- 
quency k may be used as a guide in making the selection that fits the 
trend . 


In actual computations, one further cycle of iteration beginning 
with an assumed mode given by expression (A 69 ) should be carried out to 


assess the extent to which the functions w^~"^, w^l and w^ + ^ 

a2 ’ a2 ’ ■ a2 

are free of all except the two components of the types appearing in 
equations (a6o) and (A6l). If the ratios of this cycle are not reason- 
ably constant, the unwanted components still present have to be removed 
by carrying out another cycle of iteration and again applying equa- 
tions (A64) and (a68). 


The method just described is clearly applicable in the general 


cases 





Eigenvalues having equal moduli include the special case of 
identical eigenvalues. As a basis for discussion let it be assumed that 


and that 



> c 2 = c 3 > C 4 | > . 


(A70) 


C 2 C 3 C 23 


(A71) 


The significance of the occurrence of these two identical eigenvalues 
is that the wing system may oscillate with the same frequency and arti- 
ficial damping in any of an infinite number of modes, any two of which 
axe linearly independent of each other and of the first, fourth, and 
higher modes. This infinite number of possible modes (all corresponding 
to ^ 23 ) are the infinitely many linear combinations of two basic 

linearly independent modes that are necessary and sufficient in combi- 
nation with the first, fourth, and higher modes to describe an arbitrary 
displacement of the wing system. Clearly, only two linearly independent 
modes corresponding to the double eigenvalue are required for 

analytical purposes. These two are designated w 0 and w Q as before 
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but with the reservation that. w 2 

differing linear combinations of a 
pendent modes that also correspond 


and w^ must be derivable as two 
single basic pair of linearly inde- 


Equations (A20) and (A21) are 


replaced in the present case by 


lim W a2 +1 ^ = lim C 23 Q 
n — >oo n — >oo 3 


b 2 (w 2 . wj + bj(w 3 - Wj) 


(A72) 


and 


v 


lim 

n — >oo v 


_(n+l) 

a2 


WT ' C 23 

a2 


(A73) 


Equation (A27) is replaced by 


lim w£?^ 
( n+ l) n — >oo 


lim v a2 + C 


= lim (l>2 w 2 + b 3 w 3) (A74) 


n — >oo 


'23 


The transformed second mode (equation (A72)) is in this case a 
linear combination of the first three eigenfunctions, and the so-called 
true second mode is actually a linear combination of the second and 
third eigenfunctions. 


If the iterative transformation procedure is now applied in the 
regular way to determine the transformed third mode, the third eigen- 
value, and the true third mode, the results will be as follows: The 

transformed third mode will be, like the transformed second mode, a 
linear combination of the first three eigenfunctions but will. be linearly 
independent of the transformed second mode. The so-called true third 
mode will be, like the so-called true second mode, a linear combination 
of the second and third eigenfunctions and will be linearly independent 
of the so-called -true second mode. The results will also include a 
second determination of the double eigenvalue C,^. It may therefore 

be concluded that the iterative transformation procedure is valid and 
sufficient in all cases of eigenvalue multiplicity. 
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APPENDIX B 

THE COMPLEX STIFFNESS FOR BEAMS WITH STRUCTURAL DAMPING 


The familiar concept of a complex force K(l + ig)s in simple (one- 
degree -of -freedom) vibrating systems having structural damping may be 
easily extended to continuous vibrating systems such as beams and air- 
plane wings. The quantity K is the elastic -spring constant, s is the 
displacement, Ks is the elastic -spring force, and Kgs is the 
structural -damping force. 


For a beam in flexure, the stiffness of the fibers is given by the 
modulus of elasticity E, which is analogous to the quantity K for the 
spring. The elastic stress at any point of the cross section is given 
by eE where e is the strain which is analogous to the displacement s 
Then the complex stress at any point of the cross section of a beam with 
structural damping is E(l + ig)«. The complex bending moment corre- 
sponding to this stress, obtained in the usual, way by integration of the 

moment of the stresses over the section, is El(l + ig) ^ ~ ^ 


dx' 


2 ' 


This result 


leads to the concept of a complex stiffness El(l + ig y ) for beams in 
flexural vibration with structural damping. Similarly, the complex 
stiffness of beams in torsional vibration with structural damping 
is GJ(l + ig$). The subscripts y and 0 indicate that the structural 
damping coefficient g may have a different value for torsional , 
vibrations than it has for flexural vibrations. Both g^ and g^ may 

be functions of the spanwise position x. 


1 
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APPENDIX C 

FORMULAS FOR EQUIVALENT CONCENTRATIONS 
AND INCREMENTS OF TORQUE- 

The formulas used in the numerical examples for computing equiva- 
lent concentrated loads and curvatures are those that have been derived 
in references 7 and 8. For the concentration at an end station the 
formula is 

5 1 ■ + 6p 2 ' ^ 

At an intermediate station 

S >2 = 3^(pi + 10 P 2 + ** 3 ) ( C2 ) 

1 

The significance of the quantities used in formulas (Cl) and (C2) is 
shown in the following sketch: 


(Cl) 



These formulas are based on the assumption that the distributed-load 
(or curvature) curve is a series of second-degree parabolic arcs. When 
applied to distributed flexural loads, the formulas give concentrations 
which produce the same bending moments in the wing at all the selected 
stations as the distributed load. The formulas may be correctly applied 
to distributed torsional loads only if GJ is constant over each bay. 

In this case the formulas give concentrations which produce the same 
torsional displacement at all the selected stations as the distributed 
load. For a station placed at a' discontinuity in ordinate or slope. 
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\ 


formula (Cl) must be applied to both the left and the right of the 
station and the results added. 


The formulas for obtaining increments of area beneath a curve of 
distributed torques are derived in reference 8. These formulas are 
based as before on approximating second-degree parabolas. They are , 
given here in a slightly different form which is better adapted to 
present uses. Thus 


A 1 1 + ^2 + *3) + k( q l ' 


(C3) 


*2 = !( q l + 4q 2 + q 3 ) " l( q l ‘ q 3) (C4) 

where the significance of A^^ and and of q^, q^, and q^ is 

shown by the following sketch: 



The ordinate at a discontinuity should not be used as the middle one 
of the three ordinates selected for use in formulas ( C 3 ) and (C4). The 
formulas are valid only where the three ordinates are co nn ected by a 
continuous curve . 
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ION TO OBTAIN FIRST COUPLED MODE FOR k = ». WING WITH CONCENTRATED MASS. 
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TABLE 1.- ITERATION TO OBTAIN FIRST COUPLED MODE FOR k - ». WING WITH CONCENTRATED MASS - Concluded 
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TABLE 2.- ITERATION TO OBTAIN TRANSFORMED SECOND MODE FOR k = «>. WING WITH CONCENTRATED MASS. 


fee 


Common factors for each column are given under the column headings. 


j 


3 ' — © * 1 

5 4 3 2 1 

Station 


Flexure: (a) First cycle. (b) Second cycle. (c) Third cycle. 




1 




5 

1 

2 

3 

4 

5 


1 

2 

n 


H 

6 

Station 


y (1 > 




B 

y (2) 

y (2) 


„(3) 

(3) 

y a2 


v (3) 

y(3) 

i 


1 

EC 


a2 




y a2 


1 

y a2 

; 7® 

*a2 


y a2 


1 

m 

J3T 

a2 

S? 



b 


B5P 

b 


2 

EIp. m 


B 

^0^7 2 

EIn ^ 


V)V 2 

EIu “ 

1(A) 


0 

-174.1 

174.1 

0 


0 

-4li.o 

411.0 

0 



0 

-465.8 

465.8 

0 



2 


-.281 

-107.7 

99.2 

-8.5 

30.2 

-.692 

- 252.5 

234.2 

-18.3 

26.4 


-.820 

-285.4 

265.4 

- 20.0 

24.4 


3 


-.252 

‘43.1 

33.8 

-9.3 

36.9 

-.757 

-99-7 

80.0 

-19.7 

26.0 


-.882 

- 112.0 

90.6 

- 21.4 

24.2 

24.3 

4 


-.094 

-13.1 

9.48 

-3.6 

38.3 

-.293 

- 30.2 

22.4 

-7.8 

26.6 


-.349 

-33.8 

25 . 4 

-8.4 

24.1 


5 


1 

0 

0 

0 

0 


0 

0 

0 

0 



0 

0 

0 

0 




Torsion: (a) First cycle. (b) Second cycle. (c) Third cycle. 




B 

2 

3 

4 

5 


1 

2 

3 

4 

5 


1 

2 

3 

4 

5 

6 



^a2 


^bl 

J2) 

2 

/*> 

”a2 


K 2 

J 2 > 

K 

* (2) 

^bl 

J3) 

^a2 

/ 3) 

y o2. 


,(3) 

J3) 

*b 

*u. 

„(*> 

^a2 

t w 

y B2. 

EC 

Station 


m 


7 s1 

a2 


0 (3> 

EC 



■ 

lo^7 0 
EIp 


■ 

lo S' 2 

EIn ^ 


■ 



\> k r 

Ein 

CD 2 


■ 


1.000 

13.38 

-1.10 

12.28 

12.28 


1.000 

24.94 

- 2.60 

22.34 

22.34 


1.000 


-2.95 

23.27 

23.27 




.901 

13.38 

- 1.15 

12.23 

13.58 


.997 

24.95 

- 2.72 

22.23 

22.3 


.995 

26.24 

- 3.08 

23.16 

23.3 




.624 

13-40 

- 1.28 

12.12 

19.44 


.988 

25.05 

-3.00 

22.05 

22.3 


.988 

26.40 

- 3 . ho 

23.00 

23.3 

23.3 

4 


.322 

6.70 

-.63 

6.07 

18.85 


.495 

12.51 

- 1.50 

11.01 

22.2 


.493 


-1.70 

11.50 

23.3 


5 


0 

0 

0 

0 



0 

0 

0 

• 0 



0 


0 

0 




w 


VdW\ 4 

m + £=pTj = 23 ' 8 I 5 T « 8 » 


o> 2 = 3-33 • 3 radians per second 


a2 L — . a2 j 
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TABLE 3.- ITERATION TO OBTAIN TRANSFORMED THIRD MODE FOR k = «. VINO WITH CONCENTRATED MASS. 


Common factors for each column are given under each column heading. 


Flexure: (a) First cycle. 


jl) L(l) Jl) L(l) L(2) a3 


(b) Second cycle. 


y ba2 y a3 “(IT 


EIm 


0 448.1 -44C.1 0 0 

.952 271.5 -25^.8 -11.1 5-6 
1.000 104.0 -87.0 -11.9 5.1 
.390 31.0 -24.4 -4.7 1.9 
0 0 


1 

2 

3 

4 

5 

6 

y (2) 

y a3 

(2) 

y b 

y (£) 

y bl 

y (2) 

y ba2 

y (3) 

y a3 

(3) 

~rh 

y a3 

b 


XoV 

EIu 

Cl) 2 


*O k 7 p 
EIn ^ 

0 

405.5 

-405.5 

0 

0 


1.000 

243.7 

-230.6 

-5.9 

7.2 

7.2 

.911 

91.6 

-78.8 

-6.3 

6.5 

7.2 

.339 

27.0 

-22.1 

-2.5 

2.4 

7.1 

■ 

0 

0 

am 

| 


(c) Third cycle 


v (3) v (3) (4) 

y bl y ba2 y a3 


4o^b7 2 

El u “ 


0 404.0 -404.0 0 0 

1.000 242.8 -230.0 -5.9 6.9 
.903 91.3 -78,5 -6.3 6.5 
•333 26.9 -22.0 -2.5 2.4 



Torsion: (a) First cycle. (b) Second cycle. (c) Third cycle. 
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TABLE 4.- AERODYNAMIC -INERTIA FORCE COEFFICIENTS FOR 


VARIOUS VALUES OF k FOR EXAMPLE WING 
[^Common factors for each column are given under the column headings 

Flexure 


k 

V 

P R0 

P iy 

p i0 

P Ry 

P R0 


7 2 

— 0) 

*Io> 2 

7 2 

E3 

*■07 2 

Oj 

\> b 7 p 
(XT 

0.036 

27.5 

-1444 

51.9 

-108.3 

92.5 

-75.6 

.12 

30.6 

-112.5 

13.44 

-8.27 

92.5 

-75.6 

.1443 

31.0 

- 13.2 

.10.82 

-4.14 

92.5 

-75.6 

.1590 

31.2 

-6 0.1 

9.61 

-2.53 

92.5 

-75-6 

.24 

32.0 

-23.8 

5.82 

1.35 

92.5 

-75.6 

.50' 

33.0 

-3.76 

2.39 

2.29 

92.5 

- 13.6 

OO 

33.6 

1.397 

0 

0 

92.5 

- 13.6 


Torsion 


k 

Sly 

Sl0 

Q iy 

Q I0 

Sly 

Sl0 


t>7 2 
- 6 - 0) 


t>7 2 
— 0)^ 


i 

b m 


■flfl 

. M- 

H 


■■1 

0.036 

3.67 

549 

-19.40 

68.3 

-75.6 

114.7 

.12 

2.52 

51.4 

-5.03 

11.42 

-75.6 

114.7 

.1443 

2.36 

37.5 

-4.05 

8.48 

-75.6 

114.7 

.1590 

2.28 

32.0 

-3.60 

7.24 

-75.6 

114.7 

.24 

1.98 

18.24 

-2.18 

3.67 

-75.6 

114.7 

.50 

1.623 

10.74 

-.895 

1.143 

-75-6 

114.7 

00 

1.397 

.1409 

0 

0 

-75.6 

114.7 


















TABLE 5.- ITERATION TO OBTAIN FIRST MODE FOR k = 0.1443. WING WITH CONC 








































































































































































































































TABLE 6!- ITERATION TO OBTAIN TRANSFORMED SECOND MODE FOR k = O.lMi-3. WING WITH CONCENTRATED MASS - Concluded 
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TABLE 8.- COMPUTATION OF TRUE SECOND MODE FOR k = 0.1443. 
WING WITH CONCENTRATED MASS. 

Flexural functions are in terms of bj torsional functions are in 

Cl —i 

radians. — - 1 = F 12 = 8.65 - 1.6001. 


W (4) 

* 1J \,U 


Station 

Flexure 


2 


y lN + iy ll 


F 12( y lR + y ll) 


941.0 - 165.581 | |108.8 + 0.96 




61.7 + 0.851 -8.4 + 6.31 


5.90 + 0.l8i -2.7 + 2.491 


108.8 + 0.961 
53-3 + 7.151 


21.06 + o.49i -6.9 + 6.321 14.2 + 6.811 


3.2 + 2.671 

0 


Torsion 





0W ♦ i^) . 

dIR ”bll 
F 12 (^LR + !0i^ 



2 

3 

01R + rfll 

rf (5) * -rrf (5) 
ra 2 R + i 0a2I 

-0.772 + 0.4011 

27.13 - 4.071 

-0.786 + 0.2941 

24.37 - 4.351 

-0.8l6 + 0.0481 

16.94 - 4.571 

-0.409 + 0.0261 

8.71 - 2.361 

0 

0 


02R + 1^21 






















TABLE 9.- COMPUTATION OF TRUE THIRD MODE FOR k = 0.1443. WING WITH CONC 
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Figure 1 


Transformed second mode: 

Va2 = H ' y\ 



Inertia load: 


r “ 2^( y 2 


y l) 


-ttttTTTTTTTt*. 


Forcing load: 


*{“‘2 ' “ l 2 ) S 





~ Illustration of physical basis of iterative transformation 

procedure. 
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P = 0.002378 slug/ft 3 

n = 32.6 

El = 977 .I lb-ft 2 
GJ = 48o.6 lb-ft 2 


— r— = 423,000 (radians /sec ) 2 


( EI = 

W GJ 


0.1353 


Figure 3>- Properties of cantilever wing used in numerical examples 
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Figure 


□ Experiment (reference 5) 
o Operational solution (reference 5) 
a Rayleigh-Rifz: 3 modes (reference 6) 
v Rayleigh-Ritz: 4 modes (reference 6) 
o Iteration: 4 stations 


Qa 



v, ft/sec 


(a) Variation of artificial damping with airspeed. 
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3d^ solutions 
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800 


(b) Variation of frequency with airspeed. 

4.- Variation of artificial damping and frequency with airspeed in 
first three solutions . Wing with concentrated mass . 
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Translation, y 


Rotation, <f> 





Station Station 


Figure 5.- Amplitudes and phases of modes. Wing with concentrated mass 
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□ Experiment (reference 5) 
o Operational solution (reference 5) 
a Rayleigh-Ritz: 2 modes (reference 6) 
v Rayleigh- Ritz= 3 modes (reference 6) 
o Iteration: 4 stations 


.2 
0 

9 a -.2 

-.4 
-.6 

0 200 400 600 800 

v, ft/sec 

(a) Variation of artificial damping with airspeed. 
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(b) Variation of frequency with airspeed. 

Figure 6.- Variation of artificial damping and frequency in first three 
solutions . Wing without concentrated mass . 
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